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ABSTRACT 


A computer program was developed to solve the quasi-one- 
dimensional, unsteady flow problem associated with starting 
phenomena in various nozzle-short diffuser combinations of 
interest to gas dynamic lasers. The program was based on the 
techniques developed by George Rudinger in his book, Wave Dia- 
grams for Nonsteady Flow in Ducts. The techniques were modi- 
fied to facilitate a fixed grid for computational ease. 

The program was used to analyze the starting transients 
of a Mach 3 nozzle with a semi-wedge diffuser. Results indi- 
cated that for a given inlet/exhaust pressure ratio, diffuser 
starts could be obtained at higher ramp angles for thin dif- 


fusers than for thick diffusers. 
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I. INTRODUCTION 


Gas Dynamic Lasers (GDLs) typically operate at cavity 
Mach numbers of 4.0 to 5.5 in order to obtain the kinetics 
desired for population inversion and signal gain. GDL opera- 
tion often consist of a series of short bursts which result 
in periods of supersonic flow followed by periods of no flow. 
Although efficient pressure recovery is not of prime impor- 
tance, the high Mach number, low pressure flow must be dif- 
fused in order to permit exhausting to the atmosphere. While 
rapid and efficient diffuser start-up is required in all 
laser applications due to its effect on lasing performance, 
it becomes especially critical in the short burst mode of 
operation. 

Some chemical lasers also require diffusing of low pres- 
sure high flow rate gases. Since these gases are collected 
and condensed or absorbed, diffuser efficiencies become more 
important in this type of application. 

Aircraft installation of lasers compounds the above prob- 
lems with size, weight and pressure ratio limitations. Dif- 
fusers used on GDLs currently constitute a large percentage 
of the laser size and hence become prohibitive in airborne 
installations. If an aircraft axial flow compressor is used 
for gas supply, pressure ratios are limited to approximately 
BU tot. 

The diffuser starting problem is similar to that obtained 


in starting supersonic wind tunnels. Extensive work was done 
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in this area following World War II and into the early 1950's. 
Some of the means which have been utilized in starting these 
tunnels are variable area diffusers, overpressure, perforated 
walls, auxiliary air injection, and boundary layer control 
[Ref., 6), 

In conjunction with work being done for NAVAIR at China 
Lake on the airborne gas dynamic laser, research is being 
undertaken at the Naval Postgraduate School on minimizing 
diffuser size requirements and starting time. The research 
has proceeded along “— directions; experimental and analyti- 
cal. The experimental investigation of diffuser start-up has 
been conducted by LCDR J. Zerr [Ref. 8], utilizing a blow- 
down wind tunnel to establish Mach 3 and Mach 4 flow. Asa 
first step in the analytical program it was desired to have 
a computer model by which various flow conditions could be 
analyzed. 

A complete analysis of the starting process would require 
a four-dimensional approach: 3 spatial coordinates plus the 
time coordinate. In such an approach, discontinuities and 
non-linearities caused by the presence of shock waves create 
problems of enormous complexity. 

One can obtain insight into the starting phenomenon through 
a quasi-one-dimensional analysis which utilizes one space and 
one time coordinate. Transverse area-change effects may be 
included if the radius of curvature of the change is large in 


comparison with the axial duct length. 
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Rudinger [Ref. 5] treats the one-dimensional unsteady, 
non-linear wave problem for flow in ducts. The procedures 
he develops allow for the inclusion of heat addition, vari- 
able entropy, variable area and several boundary conditions 
and discontinuities. Since his method is designed for hand 
graphical analysis, however, it does not lend itself directly 
to computer numerical methods. 

Spalding [Ref. 7] has modified the techniques of Rudinger 
in order to permit efficient digital computation procedures. 
His modifications consist primarily of the implementation of 
a fixed spatial grid combined with a backward interpolation 
of wave characteristics. Additionally, the Riemann variables 
for the characteristics have been re-defined for ease of com- 
putation. 

Solution of the characteristic-wave problem is also dis- 
cussed extenSively in Refs. 1, 3, and 6. Several alternatives 
are offered for handling the iteration process which inevit- 
ably results from the inclusion of the area effect term. The 
complexity of these schemes, in any case, must be weighed 
against the increase in accuracy and the obvious limitations 


of the one-dimensional analysis. 
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II. APPROACH 


The anita program developed in this thesis is based on 
the quasi-one-dimensional unsteady wave analysis described in 
Ref. 5 as modified by Ref. 7. These techniques have been 
utilized in the development of a previous computer program, 
PREX, for use with cyclic, constant area pressure exchangers. 
Since many of the subroutines used in PREX were general in 
nature, they were adapted as a starting point for this thesis 
and extensions were developed to permit the calculation of 
supersonic flow through ducts of varying area. 

The program essentially solves a shock-tube type problem 
in a duct with one closed end and one open end, (i.e., a 
Ludweig tube). The nozzle and diffuser are located between 
the diaphragm and the open end. When the diaphragm is broken, 
pressure waves travel through the nozzle/diffuser section 
while rarefaction waves proceed toward the closed end. Since 
the diaphragm is ddaetea closer to the open end, the diffuser 
starts prior to the time the rarefaction waves reach the 
closed end. In this manner the diffuser start-up process is 
isolated from the reflected rarefaction waves returning from 
the closed end in order to more closely approximate a reser- 
voir boundary condition. The computational model is depicted 


in Figure 1. 
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III. DERIVATION OF EQUATIONS 


In the derivations of the characteristic equations, the 

following implicit assumptions have been made: 

1. The gas flow obeys the ideal gas relations. 

2. The gas is calorically perfect (constant specific 
heats). 

3. The flow is considered quasi-one-dimensional (area 
effects are included). 

4. The duct is rigid in time with no mass flux through 


the walls (e.g., no blowing or suction). 


A. CONTINUITY 
Figure 2 depicts the mass entering and leaving a control 
volume during time dt. Conservation of mass requires: 
Increase of mass within Net mass flux through 


the control volume the control surface 


3-(pAdx)at = pAudt-[pAu + 3-(pAu)dx]dt 


3 9 » 
976 PA) + aZCpAu) = 0 


Since the duct is rigid, A is not a function of time but 
is a function only of x. Therefore expansion of the above 
equations: results in: 


oe , pudA ON gg Orie 
Ea iit c 
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Figure 2. Mass Flow Through Control Volume During Time dt. 
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Figure 3. Forces Exerted on Control Volume 
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B. MOMENTUM 
Figure 3 depicts the forces on a control surface within 
a duct of varying area. From Newton's second law: 


: _ forces acting on the body 
body acceleration “7 “mass of the body: 
Let f = summation of frictional dissipative forces/unit mass; 


then 


2 ae SEP oA sak Bis os 
Dt cade i ps5 dx [pA + ax (PA) AX] + fpAdx} 


Expanding the above equation, one gets 


du Ee eae | oA oA dp 
Go SS OA a oA op 
at * “ox ~ pAdxlPA + Pa dx - (pA + paidx + Aj ydx) + fpAdx] 
i: ee oe 
at * Uox KE. ales 
C. ENERGY 


The preceding control volume approach may also be used 
for a similar derivation of the applicable energy equation. 
Several reliable references have used this approach and the 
results will not be duplicated here; exceptionally thorough 
discussions are contained in Ref. 3 and 8. For an adiabatic, 
quasi-one-dimensional flow of a perfect gas the energy equa- 


tion reduces to: 


ce u2 ae: Me. th wcouD, . pu dA _ 
ee geey poe Yow ta pares. Ap ax 


where e is the gas internal energy/unit mass. 


D. CHARACTERISTIC EQUATIONS 
The equations of continuity, momentum and energy may be 


combined to obtain the characteristic (Riemann) wave equations. 
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Multiplication of the momentum equation by u yields 


du 200. . 3 dpe. 
= > — += = 
ot a ox p ox ou 


Subtracting the above from the energy equation, and assuming 


that f=0, gives 


The continuity equation 


So 4 pu dA, du, 99 
oo Aas ae ee 


may be re-written in the form 


9p udp, P su, pu dA _ 
ae at * ri fee Gan 


Subtracting this from the energy equation yields 


oe oe p oo udp _ 
Jk ee pa LON = OS ae 
at~ “ax "8p? at e aueee 


For an ideal gas with constant specific heats 


S) . 2 da _ dp 
d(p) y-la p 


For quasi-one-dimensional adiabatic flows the entropy of each 
particle may be expressed as 


as as 
— + — = 
t re - 
Combination with the continuity: equation produces 
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The momentum equation 


combined with the perfect gas relation 
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Addition and subtraction of Equations 1 and 2 gives 
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Note that the above two equations can be interpreted as the 
rate of change with time of the quantities (a + Yu) and 


(a - tu) along lines with slopes of +- and —— respectively 


a 
in the x, t plane. Defining new variables [Ref. 7] 
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* 
which implies 

Po=a 
the derived equations may be re-written 
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Referring to Figure 4 and applying the above equations it 


can be seen 
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which are the Riemann equations along the characteristics ex- 
pressed in terms of the newly defined variables Po: Oe» and u.. 
Reference 2 treats the constant area case of duct flow for 


which the above equations become 


Hep aa). <b a ete or, See 
pp p QQ Q pQ 2c, ox 
Bee et = Pee et fg ys SS 
SxS) 1S ss s ps a ox 


Through en analysis of wave interaction including a con- 
tact discontinuity, a justification is presented for writing 
the above equations as 


Fo Up eq" 


Sele = Poe Ue 
and o_ = that value of o determined from the particle path of 
Slope 1/u. 

Let HPS and HPQ represent the average values of 
(Sees aaa along the characteristics PS and PQ respectively. 
Then HPS dt and HPQ dt represent changes along the character- 
istics over time step dt due to the changing area. The term 


(Ses oS im may be easily obtained at points S and Q and an 


iterative procedure allows it to be evaluated at point P. 
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Figure 4, Computational Procedure. 
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The equations may now be re-written as 


P o.+U_ = P.o. +U_-HPQ dt 
poe eq q" 

Po -U 
ps p 


P_o -U _-HPS dt 
Sesoas 
which may be solved for me and U_ as 


PAd.+P_o0 +0, -U.-HPQdttHPsdt 
ss s 


ee OO Q (3) 
p OL ae 

ies Pg Pg Ug Mgt /0 , eett sat ae 
p (1/% + 1/o) 

2s = Op (5) 


where R is the x-intercept point of a characteristic having 
a Slope 1/u. 

Equations (3), (4), and (5) are the equations which are 
actually programmed. They are utilized in subroutine STEP 
and are modified by boundary conditions for use in determin- 
ing duct end values in subroutine PORTS. A discussion’ of 
these modifications as well as how the equations are utilized 


is contained in a description of the program. 
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IV. DESCRIPTION OF PROGRAM SS DIFFUSER 


A. PROGRAM EVOLUTION 

The computer program developed is based on a program 
called PREX by Lebius Matthews [Ref. 7]. PREX (mnemonic for 
pressure exchanger) solves the unsteady flow problem result- 
ing from opening and closing several valves which separate 
the duct from ports maintained at various pressure levels. 

The valves utilized are modeled as orifices which change 
in size linearly with time during the opening and closing 
process. Extensive work was done in attempting to adapt this 
concept for application in program SS DIFFUSER (mnemonic for 
super-sonic diffuser). In this application extremely rapid 
valve operation was desired in order to simulate the high 
pressure rises obtained in laser starting. Additionally, 
the valves had to be re-modeled to handle subsonic and super- 
sonic flow. 

In all our attempts to use the valves in PREX, erroneous 
stagnation pressure increases were observed near the area of 
the valve. Once established, these new stagnation pressures 
then propagated into the nozzle/diffuser section. Addition- 
ally the magnitude of the established flow Mach number in the 
duct seemed to be a function of the valve opening time. The 
valves were finally abandoned in favor of a more simplified 
boundary condition which consists of a bursting diaphragm 
such as witnessed in shock tube applications. Subroutine 


PORTS was subsequently completely re-written to handle the 
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boundary conditions as outlined by Rudinger [Ref. 5]. Minor 
modifications were made in subroutine PORTX and in the main 
program, combined with the deletion of two other subroutines 
which handled the cyclic aspect of the pressure exchanger 
problem. 

Most of the changes made to handle the variable area duct 
were made in subroutine STEP. This is also the subroutine 
which would have to be changed to include frictional or heat- 
ing effects. | 

Subroutine STEP utilizes an entropy averaging process to 
avoid the discontinuity problems associated with shock wanes, 
This is developed and discussed in depth in Ref. 6. The im- 
pact of this averaging on the results is that entropy values 
near the shock can not be relied upon to produce an accurate 
representation of the physical process. Also, if a station- 
ary shock appears somewhere in the duct, the entropy gradient 


eventually averages out entirely. 


B. EXPLANATION OF SUBROUTINES 

The following is an explanation of subroutines used in 
program SS DIFFUSER. Figure 5 contains a flow chart of the 
computational sequence. A list of symbols used in the com- 
puter program can be found in Appendix B. 

1. BLOCK DATA (data input) 

BLOCK DATA are divided into six categories of input 

data. Under FLUID PROPERTIES, the specific heat, ay and CU 


of the fluid are specified. “These are treated as constants 
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Figure 5. Flow Diagram of Program SS DIFFUSER. 
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throughout the program. Input data under BOUNDARY CONDITIONS 
are the temperatures sae pressures at the ends of the duct. 
IDATA allows one to specify static or stagnation pressure in- 
put values. In the sample program these have been non-dimen- 
Sionalized with respect to atmospheric conditions. In section 
VALVE PARAMETERS, the input value of NOVALV allows one to 
specify a fully open or fully closed right duct end. INITIAL 
CONDITIONS data fixes the values of temperature and pressure 
within the duct prior to the start of calculation. 

DUCT PHYSICAL GEOMETRY information is used by subroutine 
ARAT to determine the area effect on the flow. A more de- 
tailed discussion of these values is included in the descrip- 
tion of subroutine ARAT. PLOTTING AND PRINTING PARAMETERS 
controls the printing and plotting frequency and spacing of 
spatial nodes in the output. It also provides for a maximum 
time of solution, TMAX. 

2. MAIN (main program) 

The main program controls the order of calculations 
and specifically the calling of the designated subroutines. 
All the calculations are done on the subroutine level and 
passed to main through dummy variables and common statements. 

MAIN utilizes TMAX to determine the maximum solution 
time. 

3. INIT (initialization) 

This subroutine utilizes the specific heat inputs to 
determine the GAM- constants (see Program Terminology, Appen- 


dix B) used throughout the program. It also establishes the 
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spatial interval according to the number of calculation points 
desired. 

Slopes, areas, temperatures and pressures are initial- 
ized to the input values. From this information, initial 
values of the entropy variable, SIG, are determined. 

4. ARAT (area ratio) 

Subroutine ARAT utilizes the duct geometry to calcu- 
late the value of the area effect term in the Riemann Equa- 
tions. Its output, RAT, is a ratio of the slope of the duct 
wall to the duct area at each value of x. Areas and x-coor- 
dinates have been nCamithdens Lee iaaad with respect to duct 
length in the sample program. 

To provide for a smooth change of area in the nozzle 
section, the area change has been modeled as a sinusoidal 
variation. The sinusoid amplitude, B, is equal to one-half 
the difference between the original duct area and the area at 
the nozzle throat. Thus, a smooth contour having zero slope 
at nozzle entrance, throat, and exit is provided. 

The diffuser section allows for the modeling of a dif- 
fuser having an initial ramp input and three subsequent changes 
in flow direction. Figure 1 depicts the above models. 

5. AUX (auxiliary) 

Subroutine AUX calculates the values of output quan- 
tities other than at spatial calculation nodes. Although the 
values of A(I) and U(I) are used by subroutine STEP, other 
computations made here are not critical to calculations made 
by other routines. Their inclusion serves primarily as a 


check on input data and boundary conditions. 
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6. OUTPUT (printed output) 


This routine gathers information from other subrou- 
tines and consolidates it for printing and plotting purposes. 
The plotting routine, UTPLOT, is a catalogued routine avail- 
able under the IBM 360 system which was utilized. Similar 
such routines may be easily substituted. 

The variables shown in the sample output have been 
non-dimensionalized with respect to atmospheric conditions. 
The entropy value (s) is actually s/R. 

7. %STEP (time step) 

In a sense, STEP is the real controlling routine with- 
in the program. It controls the calculation at all interior 
points and uses routines SLOPE and LININT to do so. It es- 
sentially uses a backward linear interpolation scheme which 
may best be described by reference to Fig. 4. 

Consider the calculation at some interior point P. 

At some previous time, t,, the solution exists at points Q, 

R, and S as well as along the entire spatial grid. Since it 
is desired to maintain a constant grid interval, character- 

istics 1, m, and n are generated through P based on the pre- 
viously known values of U and A at points Q, R, and S. For 

the m characteristic shown, this would be the average of the 
values of U+A at Q and R. This same techniques is then ap- 

plied along the entire length of duct. 

The linear interpolation scheme requires that these 
characteristics intersect the x-axis at time t, somewhere be- 


tween the calculation node and the adjacent node. To 
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facilitate this, the time step is adjusted to correspond with 
the minimum slope characteristic. Additionally, since this 
characteristic is only a first approximation, only a portion 
of this step is utilized as dictated by duct geometry and the 
expectation of area effect. It may be readily seen that the 
calculation time is quite sensitive to the number of grid 
points. 

Once the time interval has been established, the cor- 
responding x-intercepts are calculated by subroutine SLOPE. 
Values of the computing variables are then calculated at these 
intercepts by linear interpolation of the bracketing calcula- 
tion nodes Q, R, and S. Using these values, a rough guess 
for the calculation variables are obtained for point P based 
on the Riemann equations with no area change. 

The area effect term is integrated along the charac- 
teristic by assuming it to be a constant equal to the average 
of the end values. Subsequently new values are obtained at 
P through the introduction of this term. An iteration pro- 
cedure is introduced here in order to obtain a more accurate 
area ,cifect, term. p 

Upon completion of the velocity and pressure calcula- 
tion at point P, the entropy term, SIG, is updated at time 
Bit, dt. 

8. SLOPE (characteristic slopes 

As previously mentioned, subroutine SLOPE is respon- 

sible for computing the slopes of the characteristics. Char- 


acteristics which pass through endpoints of the duct are 
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calculated only for the one originating within the duct. 
These are then combined with boundary values in subroutine 
PORTS (Fig. 4). 

Two options may be used in calling SLOPE. The first 
option computes a slope on the basis of the average values 
of the bracketing spatial nodes. This is the option which is 
utilized in the sample program. The other option, which may 
be used with appropriate changes in subroutine STEP, allows 
the values in the slope equation to be externally dictated. 

In either of the above options, the subsequent x- 
intercept at time t-dt is computed based on the slope and 
time interval. 

9. LININT (linear interpolation) 

LININT performs a simple linear interpolation for 
variables based on the values of the x-calculation nodes 
bracketing the desired spatial coordinate. The interpolation 
always takes place between two adjacent nodes. Endpoint ex- 
clusion is accomplished in the Same manner as described for 
subroutine SLOPE. 

10. PORTX (port auxiliary) 

PORTX is the controlling routine for calculation of 
pressure and temperature variables external to the duct. If 
the duct ends are closed, the dummy variables passed to PORTS 
are appropriately set to zero. If the duct end is open, as 
in the case of the sample program, the variables are calcu- 
lated based on whether static or stagnation values have been 
specified. These values are then passed to PORTS where they 
are meshed with the characteristic values. 
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11. PORTS (port solutions) 


Subroutine PORTS controls the calculation of the ve- 
locity and pressure variables at the duct ends. If the ends 
are closed, PORTS provides the appropriate wave reflection 
from the Riemann equations. 

For outflow, PORTS differentiates between two cases, 
sub-critical (subsonic) and super-critical (supersonic). For 
the sub-critical case the existing outside pressure is uti- 
lized as discussed in Ref. 5. For supersonic flow, all char- 
acteristics dkiaiiase from within the duct and hence the 
Riemann equations are directly applicable. PORTS also pro- 
vides an indicator to OUTPUT indicating whether the flow is 


supercritical or subcritical. 
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V. DISCUSSION OF RESULTS 


The computer program was initially used to solve a sample 
problem from Ref. 5. In Appendix A the computed solution is 
presented along with the graphical solution given by Ref. 5. 
The results can be seen to be in very good agreement with 
those of Rudinger. 

Additional flows were considered utilizing nozzle con- 
traction ratios of 4:1 and gas specific heat values Cy and c 
of 0.172 and 0.24 BTU/lbm-°R, respectively. Isentropic values 
for these conditions indicate a resulting Mach number of 2.94 
should be obtained. Mach numbers of approximately 3.04 were 
obtained indicating an error of 3.4%. This resulting error 
is quite low considering the inherent limitations of the lin- 
ear interpolation process utilized. 

Tables I-V depict typical results obtained for a driver 


pressure value below that required to obtain a fully started 
a,t 
lig 
foot duct length. The pressure utilized was p/p) of 16.0 


condition. Time was nondimensionalized ( ) utilizing a four 
(po = ambient pressure). A semi-wedge diffuser (see Fig. 6) 
was utilized having a ramp slope of 0.05 (6=3°) and contrac- 


ee 
_aitiuser throat) or"0. 94. 


tion ratio ( Subsonic flow is es- 


test section 


tablished ahead of the nozzle throat (x=.8125) with supersonic 
flow in the diverging section. A formation of a stationary 
shock wave (Fig. 6) can be observed in the diverging section 
(x=.8700) of the nozzle resulting in subsonic flow in the 


test section (x=.8750 to x=.8900). 
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Tables VI-IX depict values for the same duct geometry for 
which the driver pressure has been elevated to 20.0 to drive 
the shock wave out of the nozzle and downstream (Fig. 7). 

Pressures necessary for starting were considerably greater 
than witnessed in the laboratory. Possible explanations for 
this are that the program incorporates no friction, the burst- 
ing diaphragm causes a stagnation pressure loss, and only nor- 
mal shocks are modeled. An effective stagnation pressure 
ratio of approximately 18.0 can be observed in the nozzle. 
When combined with the losses across a Mach 3.0 shock 
(Pyo/P447-33), a correlation can be seen with laboratory ob- 
served starting values of 6.0. 

These results appear to be in good agreement with the 
cited references. References 1 and 4 discuss similar condi- 
tions extensively. Although they treat the diverging section 
as a step increase in duct area, antidise results are reported. 
For weak shock waves passing through a widening duct, a re- 
flected expansion wave is created along with a transmitted 
shock. Additionally, a contact discontinuity is created by 
the entropy increase due to the shock. The resulting flow is 
subsonic in all regions; the entire process is shown in Fig- 
ure 8A. For stronger shock waves region 5 disappears (Figure 
8B), and at yet stronger values a stationary shock occurs as 
shown in Figure 8C. Further increases in arriving shock 
strength cause the stationary wave to be swept downstream as 
depicted in Figure 8D. Although these distinct waves re 
be easily picked out in the computer results due to complica- 
tions caused by a continuously varying area and the effects 
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of the nozzle contraction preceding in the area of interest, 

the general trends described above are in qualitative agree- 
ment. It appears as though the area contraction has a strength- 
ening effect on the shock and causes an increase in shock | 
velocity as might be expected from continuity. This even- 
tually results in the formation of a reflected nave travel- 

ling upstream. Such an effect can be seen when a diffuser is 
added to a nozzle with the resulting larger pressure gradient 
across the stationary shock wave in the nozzle diverging sec- 
tion. 

Figure 9 shows the effect of diffuser contraction ratio 
and ramp angle on start-up. For a fixed nozzle entrance to 
diffuser exit pressure ratio of 20.0, startups were attempted 
utilizing 56 combinations of ramp angle and contraction ratio. 
The results indicate that thin diffusers permit the use of 
higher ramp angles than do thick diffusers. Zerr [Ref. 8] 
determined similar results experimentally and advocates the 
use of thin, high ramp angle diffusers for efficient starting. 
Values for high ramp angles and extremely thin diffusers are 
not shown since these values slip through the computational 


grid size utilized (x=.01). 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The developed computer model appears to give satisfactory 
results in one-dimensional wave analysis of the diffuser 
starting phenomena. Only a few of the possible parameters 
have been varied thus far but the program may be used for 
more extensive studies in.-a-number of areas. 

The high contraction ratio nozzle did not appear to vio- 
late the implicit quasi-one-dimensional assumption. Reference 
1 gives the criterion that the axial component of the flow 
velocity and its derivatives should be larger than the cor- 
responding transverse component by at least one order of mag- 
nitude, to satisfy this assumption. While this criterion may 
be overly restrictive in some applications, it must be abided 
by in high contraction ratio problems. 

The program appears to have no difficulties in the tran- 
sonic range where one of the characteristics develops a ver- 
tical slope. A problem was anticipated in this regime, and 
slope and interpolation equations were written so as to avoid 
indeterminate forms. 

One of the most restrictive features of the computer pro- 
gram is the time step dependence on grid size. The interpo- 
lation process requires that the characteristics at time step 
t + dt originate within one spatial grid space at time t. 
Thus in decreasing the grid size for more accuracy one must 


pay a high price in calculating time. 
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As previously mentioned, extensions may be made in sub- 
routine STEP to include frictional and heat transfer effects. 
Additional changes may be possible to include other boundary 
layer phenomena. 

Studies so far have been based on diffusers with a lead- 
ing ramp angle followed by a constant area section to the duct 
end. A similar study may be done on a full wedge diffuser 
(converging/diverging) in order to gain a more complete pic- 
ture of diffuser effects. 

Reference 2 documents the results of several applications 
of the basic techniques utilized by the program. These in- 
clude a shock tube, supercharger and pressure exchanger. 
Program SS DIFFUSER may be used for constant area shock tube 
problems by setting the diffuser wedge angles and nozzle 
sinusoid equal to zero and utilizing the closed end condition 


for the right end. 
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APPENDIX A 


RUDINGER DUCT PROBLEM 


The following problem was taken from Ref. 4. The compu- 
ted output is presented for comparison with the graphical 


solution presented by Ref. 4. 


PROBLEM : 

A duct of constant cross section is filled with air 
(y=1.4) that is isentropically compressed from atmospheric 
pressure to P,=2.83. The duct is suddenly opened to the at- 
mosphere. How does the pressure vary at the closed end of 


the duct? 


SOLUTION: 

The pressure history at the duct end is depicted in Fig- 
ure Al. Close agreement is seen throughout. Initially an 
expansion wave travels toward the closed end and is reflected 
back to the open end. The open end condition causes it to 


reflect back as a shock wave which reaches the closed end 


shortly after a = 30, 

There is a slight variation of results in the vicinity of 
the shock wave where the computer solution depicts a finite 
slope instead of the predicted vertical discontinuity. This 


is a result of the averaging process across the shock previ- 


ously mentioned. 
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Fortran Symbol 
A(I) 
AA 
AK 
AM(I) 
AMAC(I) 
AM1 


AM1(I) 


AM2(I) 


AMT 


AN(I) 
AN1 


AN1II 


AN1(1) 


AN2(T) 


AREA 


ANT 


APPENDIX B 


PROGRAM TERMINOLOGY 


Meaning 


Sonic velocity 

Duct area 

Fraction of maximum allowable time step 
Velocity variable 

Local Mach number=U/A 

Initial 'AM' in duct 


'"AM' interpolated along 'N' character- 
istic 


'AM' interpolated along 'M' character- 
istic 


'AM' prior to inciusion of area effect 
term 


Pressure variable 


Initial value of 'AN' in duct 
(low pressure side) 


Initial value of 'AN' in duct 
(high pressure side) 


"AN' interpolated along 'N' character- 
istic 

"AN' interpolated along 'M' character- 
istic 

Duct area (in ARAT) 


'AN' prior to inclusion of area effect 
term 


Subroutine name 


Amplitude of nozzle sinusoidal area 
variation 
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Bt 
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BE(I) 
BEM(I) 
BEN(1) 
BES(I) 
‘ BMAX 
Cc 

C1 

ce 

C3 
CELL 
COMP 
cP 

cv 

CYC 


DADX 


DG 
DL 
DT 
DX 
EPS 
FLOW 
FRA 


G(NPORT) 


GAM 


GAMC 


Characteristic slope 

Slope of 'M' characteristic 
Slope of 'N' characteristic 
Slope of 'S' characteristic 
Maximum of BEM, BEN, and BES 
Label for common block 

Label for common block 


Label for common block 


‘Label for common block 


Label for common block 

Label for common block 

Specific heat at constant pressure 
Specific heat at constant volume 
Label for common block 


Da/dx change in duct area with respect 
to x 


Increment in 'G' 

Duct length 

Time step size 

Spatial step size 

Iteration control term 

Label for common block 

Fraction in linear interpolation 


Integrated mass flow into port numbered 
'NPORT' 


Ratio of specific heats=CP/CV 


Critical pressure ratio for 
INFLOW=SQRT( (GAS+1)/2) 
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(TAOTH)D Te 


Mie 
OAD 


HPQ 


IDATA(I) 


IHAND 


IND 


INDL 
INDR | 
INIT 


IOPEN 


K1 


INDEX 


Label for common block 

= (GAM-1)/2 

= -2/GAM1 

= 2(GAM-1)/(GAM+1) 

= (2/(GAM+1))**((GAM+1) /(2*(GAM-1))) 
= GAM/GAM1 

= 1/GAM5 

= (GAM+1)/(2+GAM) 

= (GAM+1)/(GAM-1) 


Area effect term along 'M' character- 
istic 


Area effect term along 'N' character- 
istic 


1 If port static pressure supplied 


2 If port stagnation pressure supplied 


Index identifying port 1 for left hand 
port, and 2 for right 


Indicator to route computation 1 for NC 
flow----- closed port 8 for super-criti- 
cal outflow through fully open port 9 
for sub-critical outflow through fully 
open port 

‘IND' at left end 

'IND' at right end 

Subroutine name 

Index, =2 if end closed, =1 if end open 
Alphameric symbol identifying end 
Alphameric symbol for right end 

Counter 


Counter in DO loop 


KDPM + 1 
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Ane FEN) y it ps 


KDPM Location of diaphragm in duct 


L Causes print of variables every L-TH 
node 

LININT : Subroutine name 

MAIN The main program 

N Numbers of grid intervals in spatial 
direction 

N1 Numbers of grid points in spatial 


direction = Nt1 


NOPENL Number of port open on left hand side 

NOPENR Number of port open on right hand side 

NOVALV NOVALV = 1 Implies right end open, 
closed for all other values 

NPLOT Plot control 

NPRINT Print control 

NSTEP Number of time steps 

OUTPUT Subroutine name 

P(I) p Pressure 

PHI Duct end opening (if open, if closed) 

PHI Sinusoidal angle in ARAT 

PHIL '"PHI' at left port 

PHIR 'PHI' at right port 

PI Initial pressure in duct (low pressure 
side) 

PII Initial pressure in duct (high pressure 
side) 

PIE 3.1415926536----- 

PO Pot Stagnation pressure 

POCL Stagnation pressure in duct, left 


boundary 
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POCR 


POL 


POP(I) 


POR 


PORTS 
PORTX 
POS 
PR 
PRL 
PRR 
PROP 
PX 
PXL 
PXP(I) 
PXR 
Pxs 
RAN(I) 
RANGE 
RATIO 


RAT(I) 


RATX(1) 
RATX1(1) 
RATX2(I) 
RO(I) 


eo) 


ex 


Stagnation pressure at duct, right 
boundary 


Stagnation pressure at left end (exter- 
nal) 


Stagnation pressure at duct end (exter- 
nal) 


Stagnation pressure at right end 
(external) 


Subroutine name 

Subroutine name 

Non-dimensional spatial coordinate 
Ratio of stagnation to static pressure 
"PR' at left end (external) 

'PR' at right end (external) 

Label for common block 

Static pressure (external) 

Static pressure at left end (external) 
Static pressure in port 'I' 

Static pressure at right end (external) 
Dummy argument for 'PXL' or 'PXR' 
Limits of output plot 

Label for common block 

Label for common block 


Dummy variable in ARAT for 


RAT evaluated at C(I) 


RAT interpolated along 'N' characteristic 


RAT interpolated along 'M' characteristic 


Specific density 


Entropy 
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‘sQ) 
S1(I) 
SIG(I) 
SIG1(I) 
SIG2(I) 
SIGI 
SIGII 
SL 
SLOPE 
SOL 
SOR 
START 
STEP 
T(I) 
THETA(I) 


foe 
fit 


TIM 
TMAX 
TO 


TOCL 
TOCR 
TOL 


TOP(I) 


Dummy variable in linear interpolation 
Dummy variable in linear interpolation 
Dependent variable, function of entropy 
SIG interpolated along 'N' characteristic 
SIG interpolated along 'M' characteristic 
Initial value of SIG (low pressure side) 
Initial value of SIG (high pressure side) 
Slope in linear interpolation 

Subroutine name 

Stagnation entropy at left end (external) 
Stagnation entropy at right port 

Label for common block 

Subroutine name 

Temperature 

Angles in diffuser geometry 


Initial temperature in duct (low pres- 
sure side) 


Initial temperature in duct (high pres- 
sure side) 


Normalized time 
Cycle time 
Stagnation temperature 


Stagnation temperature in duct, left 
boundary : 


Stagnation temperature in duct, right 
boundary 


Stagnation temperature at left end 
(external) 


Stagnation temperature at duct end 
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Stagnation temperature at right end 
(external ) 


Velocity 

Dummy variable 

Spatial distance 

X interpolated along BEN 
X interpolated along BEM 


X interpolated along BES 
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